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^ ; ABSTRACT 

We demonstrate the use of the 3D Monte Carlo radiative transfer code phaethon to 
model infrared-dark clouds (IRDCs) that are externally illuminated by the interstellar 
radiation field (ISRF). These clouds are believed to be the earliest observed phase of 
Q ■ high-mass star formation, and may be the high-mass equivalent of lower-mass prestel- 

^ ' lar cores. We model three different cases as examples of the use of the code, in which 

C/^ . we vary the mass, density, radius, morphology and internal velocity field of the IRDC. 

, ■ We show the predicted output of the models at different wavelengths chosen to match 

the observing wavebands of Herschel and Spitzer. For the wavebands of the long- 
wavelength SPIRE photometer on Herschel^ we also pass the model output through 
^ ■ the SPIRE simulator to generate output images that are as close as possible to the 

I ones that would be seen using SPIRE. We then analyse the images as if they were real 

On ■ observations, and compare the results of this analysis with the results of the radiative 

CO ■ transfer models. We find that detailed radiative transfer modelling is necessary to ac- 

curately determine the physical parameters of IRDCs (e.g. dust temperature, density 
\^ ■ profile). This method is applied to study G29. 55+00. 18, an IRDC observed by the 

' Herschel Infrared Galactic Plane survey (Hi-GAL), and in the future it will be used 

. to model a larger sample of IRDCs from the same survey. 
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1 INTRODUCTION analogues to prestellar cores (Ward- Thompson et al. 1994, 

^ ^ 1 i-rr- 1 r 1 2007), fiom wMch lower-mass stars form (e.g. Solar- type 

High-mass stars evolve somewhat dirterently from low-mass ^ \ t r x tt^t^^ i xi ^ n 

_ 1 / 1 \ r . 1 . stars), m tact, IKDCs are so dense that they are optically 

stars. There appear to be (at least) tour mam evolutionary ^i • i ^ • r i i ^i / o 0.1 \ t r 

r , . , r • / ry. 1 o ^ thicK cvcu at lutra-rcd wavelengths (^8-24/imj. m many ot 

staeres ot high-mass star tormation (e.g. Zmnecker & Yorke ^1 1 1 ^1 • 1 ^ ^ i 1 ^ r 

° . -r iiii 11 11 these clouds there is little or no detectable star tormation 
2007): mtra-red-dark clouds; hot molecular cores: compact / t-» ^ 1 or^r^r^\ 
TTTT . /. 1 1. 1 11 TTTT (^-g- Parsous et al. 2009). 
Mil regions (including hyper- and ultra-compact Mil re- 
gions); and classical HII regions. The mass of an IRDC can be hundreds, or even thou- 
The earliest observed phase of high-mass star formation sands of solar masses (e.g. Rathborne et al. 2009). Often, 
is the infra-red-dark cloud (IRDC) stage. IRDCs are very dense, dark cores can be seen within each cloud. These are 
dense, massive interstellar clouds in which high-mass stars known as infrared-dark cores, and can have masses of up 
are believed to form. They were first identified in infra-red to ^100 M©, and are typically only 0.1 pc or less in ra- 
surveys, by the Infrared Space Observatory (ISO; Perault dius. Hence they have densities of up to ^10^^ hydrogen 
et al. 1996) and the Midcourse Space Experiment (MSX; molecules m~^ (e.g. Egan et al. 1998; Carey et al. 1998, 
Egan et al. 1998), as denser, higher-mass and more distant 2000). These are the most likely candidate sites of high-mass 

star formation. 

- E-mail:D.Stamatellos@astro.cf.ac.uk The Herschel Space Observatory was launched on 2009 

t Herschel is an ESA space observatory with science instruments May 14, and carries a 3.5-metre diameter telescope, and 

provided by European-led Principal Investigator consortia and three scientific instruments designed to carry out imaging 

with important participation from NASA. and spectroscopy (Pilbratt et al. 2008, 2010). The Spectral 
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and Photometric Imaging Receiver (spireEI) is the long- 
wavelength camera on board Herschel (Griffin et al. 2009, 
2010). The Hi-GAL survey is a Herschel Open Time key 
project that is mapping the inner Galaxy in five bands be- 
tween 60 and 500 /xm (Molinari et al. 2010a,b). Hi-GAL 
will map the area common to the Spitzer GLIMPSE & 
MIPSGAL (Rieke et al. 2004) and MSX (Price et al. 2001) 
Galactic plane surveys, but at FIR and submm wavelengths. 
Hence, the IRDCs will be seen in emission, in contrast with 
the previous surveys. One of the science goals of Hi-GAL 
is to determine the physical properties of IRDCs and hence 
constrain the initial conditions of high-mass star formation. 

In this paper we demonstrate the use of the 3D Monte 
Carlo radiative transfer code phaethon in converting two- 
dimensional imaging data of IRDCs from the Hi-GAL sur- 
vey, into the physical parameters of the clouds being stud- 
ied, with particular emphasis on the data from SPIRE. In 
Section 2 we describe the radiative transfer code and the 
specific input parameters used in this study. In Section 3 we 
present radiative transfer models of three IRDCS with dif- 
ferent morphologies (spherical, flattened, and turbulent). In 
Section 4 we use the SPIRE photometer simulator to pro- 
duce simulated observations of the IRDC models, and in 
Section 5 we compare the radiative transfer models with the 
simulated observations. In Section 6, we apply this method 
to determine the properties of G29. 55+00. 18, an IRDC from 
the Hi-GAL survey. Finally, in Section 7 we summarize our 
results. 



2 NUMERICAL METHOD: MONTE CARLO 
RADIATIVE TRANFER 

2.1 Code overview 

The radiative transfer calculations are performed using 
PHAETHOfJfl, a 3D Monte Carlo radiative transfer code (Sta- 
matellos 2003; Stamatellos & Whitworth 2003, 2005; Sta- 
matellos et al. 2004, 2005, 2007). The code uses a large num- 
ber of monochromatic luminosity packets to represent the 
radiation sources in the system (stars and/or the ambient 
radiation field). The luminosity packets are injected into the 
system and interact (are absorbed, re-emitted, scattered) 
stochastically with it. If a luminosity packet is absorbed its 
energy is added to the local region and raises the local tem- 
perature. To ensure radiative equilibrium this packet is re- 
emitted immediately with a new frequency chosen from the 
difference between the local cell emissivity before and after 
the absorption of the packet (Bjorkman & Wood 2001; Baes 
et al. 2005). 

The input parameters of the code are: (i) the density 

^ SPIRE has been developed by a consortium of institutes 
led by Cardiff University (UK) and including Univ. Lethbridge 
(Canada); NAOC (China); CEA, LAM (France); IFSI, Univ. 
Padua (Italy); lAC (Spain); Stockholm Observatory (Sweden); 
Imperial College London, RAL, UCL-MSSL, UKATC, Univ. Sus- 
sex (UK); and Caltech, JPL, NHSC, Univ. Colorado (USA). This 
development has been supported by national funding agencies: 
CSA (Canada); NAOC (China); CEA, ONES, CNRS (France); 
ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC (UK); and 
NASA (USA). 

^ http:/ /www. astro. cf.ac.uk/pub/Dimitrios.Stamatellos/phaethon 



Table 1. IRDC parameters used in the radiative transfer models. 
Uc- cloud central density; ro: cloud flattening radius; i^cloud- cloud 
extent; Mdoud: cloud mass; ddoud- cloud distance. 





Spherical 


Flattened 


Turbulent 


He (cm~3) 


5 X 10^ 


10^ 


5 X 10^ 


ro (pc) 


0.08 


0.04 


0.08 


^cloud (pc) 


1 


1 


1.5 


Meloud (Mo) 


1300 


510 


27200 


c^cioud (kpc) 


2 


2 


2 



profile of the system, (ii) the dust properties (e.g. opacity), 
and (iii) the radiation sources present (e.g. ambient radiation 
field and/or stars). The code calculates the dust tempera- 
ture, and produces observables: (i) the SED of the system at 
different viewing angles, and (ii) synthetic images of the sys- 
tem at different viewing angles and at different wavelengths. 



2.2 External radiation field and dust properties 

A starless IRDC is assumed, hence the cloud is only heated 
by the ambient radiation field. For the ambient radiation 
field we adopt a revised version of the Black (1994) inter- 
stellar radiation field. This consists of radiation from giant 
stars and dwarfs, thermal emission from dust grains, cos- 
mic background radiation, and mid-infrared emission from 
transiently heated small PAH grains (Andre et al. 2003). 
This radiation field represents well the interstellar radiation 
field in the solar neighbourhood but it is probably an un- 
derestimate of the ambient radiation field in the Galactic 
Plane. An enhanced ambient radiation field could result in 
higher dust temperatures at the surface of the clouds but 
the temperatures in the central cloud regions are not ex- 
pected to be significantly higher (only a few degrees). The 
luminosity packets representing the ambient radiation field 
(typically a few 10^° packets) are injected from the outside 
of the cloud with injection points and injection directions 
chosen to mimic an isotropic radiation field incident on the 
cloud (Stamatellos et al. 2004, 2007). 

The dust grains in dense clouds are expected to coagu- 
late and accrete ice mantles, thus we use the Ossenkopf & 
Henning (1994) opacities for the standard Mathis, Rumpl 
& Nordsieck (1977) interstellar grain mixture (53% silicate 
& 47% graphite), with grains that have coagulated and ac- 
creted thin ice mantles over a period of 10^ years at densities 
lO^cm"^. We assume a gas-to-dust mass ratio of 100. 



3 RADIATIVE TRANSFER MODELS OF 
IRDCS 

We model 3 different cloud geometries: (i) a spherical cloud; 
(ii) a flattened cloud; and (iii) a turbulent cloud with sub- 
structure. We also vary the mass, density and radius. The 
parameter values are chosen so as to illustrate some of the 
variety of parameter values that the code can handle. The 
parameters of each of the three cases are listed in Table 1. 
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Figure 1. Dust temperature profile of a spherical IRDC. 
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Figure 2. SED of a spherical cloud at a distance of 2 kpc. 



3.1 Spherical IRDC 

For the spherical cloud we adopt a Plummer-like density 
profile (Plummer 1915), 



n(r) 



(1) 



where ric is the density at the centre of the cloud, and ro 
is the extent of the region in which the density is approxi- 
mately uniform. This density profile is consistent with ob- 
servations of low-mass cores in nearby star forming clouds, 
such as L1544, L63 and L43 (e.g. Kirk et al. 2005, 2006, 
2007), and has been used to model such clouds (Whitworth 
& Ward-Thompson 2001). The values of these parameters 
are given in Table 1. The mass of the cloud is 1300 Mq. 

Figure 1 shows that the temperature in the cloud drops 
from 16 K at the edge of the cloud to 4 K at the centre 
of the cloud. The central dust temperature is 2 — 3 K lower 
than the typical temperatures in the central regions of low- 
mass (i.e. a few solar masses) clouds (~ 7 K; e.g. Stamatellos 
& Whitworth 2003; Stamatellos et al. 2004). This is due to 
the much higher mass and extinction of the IRDC. Heating 
from cosmic rays, which is not included in our models, may 
increase the temperature in the centre of the cloud to ~ 
5 K. Additionally, a stronger external radiation field in the 



20 



^ 15 



10 



0) 



03 



0.01 



0.1 
r (pc) 



Figure 3. Dust temperature profile of a flattened cloud at two 
different directions in the cloud. The bottom line corresponds to 
the midplane of the flattened structure, whereas the upper line 
corresponds to the direction perpendicular to the midplane. 
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Figure 4. SED of a flattened cloud at a distance of 2 kpc. 



Galactic plane could also possibly increase the temperature 
by a few degrees). 

Figure 2 shows that the SED of the cloud peaks at 
around 200 ^im, consistent with the low temperatures in 
the cloud. Because the temperature gradient in the IRDC 
is larger than the temperature gradient in a low-mass core, 
the SED is broader. 



3.2 Flattened IRDC 

For the flattened cloud we use a density profile of the form 

\ 2 



n(r, 0) = ric 



l + ^(^) sir^^ie) 



(2) 



where ric is the density at the centre of the cloud, and ro 
is the extent of the region in which the density is approx- 
imately uniform (see Stamatellos et al. 2004). The values 
of these parameters are given in Table 1. The parameter A 
determines the equatorial-to-polar optical depth ratio e , i.e. 
the maximum optical depth from the centre to the surface 
of the cloud (which occurs at ^ = 90°), divided by the min- 
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Figure 5. Density on the y = —0.2 pc (left), y = pc (centre) and y = +0.2 pc (right) planes for a turbulent IRDC. The density is 
given in units of g cm~^. 



imum optical depth from the centre to the surface of the 
cloud (which occurs at ^ = 0° and = 180°). The param- 
eter p determines how rapidly the optical depth from the 
centre to the surface rises with increasing 0, i.e. going from 
the north pole at ^ = 0° to the equator at ^ = 90°. In this 
model we assume e = 1.03, i.e. a slightly flattened cloud, 
and p — 4. This geometry may be more realistic than the 
spherical cloud presented in the previous section. The mass 
of the cloud is 510 Mq (see Table 1). 

Figure 3 shows that the temperature profile in the cloud 
is similar to the case of the spherical cloud; the temperature 
drops from 18 K at the edge to 5 K in the centre of the 
cloud. However, due to the flattened cloud geometry the 
cloud 'equator' is colder than the cloud 'poles'. 

Figure 4 shows that the SED of the cloud is also similar 
to that of the spherical cloud. There is no dependence on 
the viewing angle despite the fact that the optical depth to 
the centre of the cloud becomes Tdoud 1 at ^ 200 ^im. 
This is because the temperature varies with the direction 
in the cloud only in the outer region of the cloud, which is 
optically thin even at wavelengths up to ^ 200 ^im. 



3.3 Turbulent IRDC 

We finally examine a more realistic cloud geometry, i.e. a 
turbulent cloud. To produce this cloud we start off with a 
spherical cloud having a density profile 

n(r) = Tic , (3) 

where ric is the density at the centre of the cloud, and ro 
is the extent of the region in which the density is approxi- 
mately uniform. Note that this profile is less steep than the 
profile used for the spherical IRDC, hence the higher mass 
of this cloud. 

The values of the parameters are again given in Ta- 
ble 1. The mass of the cloud is 27200 Mq. A large turbu- 
lent velocity field (a = 0.5) is imposed to the cloud (e.g. 
Goodwin et al. 2004), and the cloud is evolved using the 
SPH code DRAGON, for enough time to produce substruc- 
ture, i.e. until cores form in the cloud (e.g. Stamatellos et 
al. 2007). DRAGON is a gravitational hydrodynamics code 
which invokes a large number of particles to represent a 
physical system. It uses an octal tree (to compute gravity 
and find neighbours), adaptive smoothing lengths, multiple 



particle timesteps, and a second-order Runge-Kutta inte- 
gration scheme. The resulting cloud is shown in Figure 5, 
in which we plot its density on 3 planes (y — —0.2 pc, left, 
2/ = pc, centre and y = +0.2 pc, right). 

We perform a radiative transfer simulation for this 
cloud using the method of Stamatellos & Whitworth (2005). 
The calculated dust temperature is presented in Figure 6. 
The dust temperature is similar to the temperature calcu- 
lated in the previous cases. However, the temperature distri- 
bution at a particular distance from the centre of the cloud 
is broader due to the dumpiness of the cloud. Figure 7 shows 
that the SED of the cloud peaks at longer wavelengths than 
in the previous two cases, as the cloud is more massive and 
consequently it is cooler. 



4 SPIRE SIMULATED OBSERVATIONS 

Simulated observations of the modelled IRDCs are per- 
formed, to assess the impact of instrument systematics. The 
simulations are produced using V2.31b of the SPIRE pho- 
tometer simulator (Sibthorpe, Chanial and Griffin 2009). 
This software simulates both the SPIRE instrument, and the 
Herschel observing modes. The simulations are performed in 
a configuration which matches the Hi- GAL observations. 

The fast scanning speed (GO^^per second) option is used 
(see the SPIRE Observers' Manuajfl), and the observations 
are made in the SPIRE parallel mode, to match the Hi- 
GAL mapping mode. These observations are performed by 
continuously scanning the telescope backwards and forwards 
in a raster pattern while the instruments continuously take 
data. 

For Hi- GAL two maps are obtained by scanning first 
along one axis of the instrument array, and then scanning 
along the perpendicular axis of the array. This provides re- 
dundancy in the data, allowing for the use of maximum like- 
lihood techniques to be used in the subsequent data reduc- 
tion (Sibthorpe et al. 2008). 

Realistic noise with a 1// spectrum is simulated with 
each of the SPIRE bolometers using independent noise pa- 
rameters based on pre-flight test data. Thermal drifts are 
not included in these simulations as the information is 



^ http:/ /herschel. esac.esa.int/Documentation.shtml 
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Figure 6. Dust temperature of a turbulent IRDC. The points 
correspond to different locations in the cloud. 



10-« F 




A(^m) 

Figure 7. SED of a turbulent IRDC at a distance of 2 kpc. 

not available to model these accurately. However, the Her- 
schel/SFlRE pipeline will automatically remove these drifts, 
making their omission from these simulations insignificant. 

The simulated data are calibrated and reduced from 
time-line data to maps using a simple map-maker. Maps at 
the three SPIRE bands, 250, 350 and 500 /xm are obtained, 
having full width half maximum beams sizes of 18, 25 and 
36 arcseconds respectively. 

5 COMPARISON OF IRDC MODELS TO THE 
SIMULATED OBSERVATIONS 

5.1 Maps 

Synthetic maps of the plane-of-sky intensity distribution of 
the three IRDCs are produced using the radiative transfer 
model. The maps are produced at eight wavebands (8, 15, 
70, 160, 250, 350, 500, and 500 /xm) that replicate the filter 
sets of recent telescope facilities. The 8 and 15 /xm wave- 
lengths match two bands from the MSX satellite (Price et 
al. 2001). The near infra-red 8 /im is also found on Spitzer's 
IRAC camera (Fazio et al. 2004). The 8 /xm filter can show 
strong emission from PAH grains (Flagey et al. 2006). The 
radiative transfer model does not include transiently heated 



small grains. These grains dominate the emission at some 
IR bands (e.g. Li & Draine 2003). In this paper we focus at 
longer wavebands (A > 70/im), where there is no emission 
from such grains. 

The far- infrared bands at 70 and 160 /im bands match 
approximately the 70 and 160 /xm bands of Spitzer's MIPS 
camera (Rieke et al. 2004), the N60 and N160 bands of 
AKARI (Kawada et al. 2007), and the 70 and 160 bands 
of HerscheVs PAC^I camera (Poglitsch et al. 2009, 2010). 
The submillimetre bands at 250, 350, and 500 /im match 
approximately the filter set of HerscheVs SPIRE (Griffin et 
al. 2009, 2010) camera and the BLAST balloon experiment 
(Pascale et al. 2008). The 850 /xm band matches approxi- 
mately the long-wavelength filter of SCUBA (Holland et al. 
1999) and SCUBA2 (Holland et al. 2006) on the JCMT. 
The 500 /im band is also roughly equivalent to SCUBA's 
and SCUBA2's 450 /xm band (Holland et al. 1999; 2006). 

Given a synthetic map of intensity it is possible to sim- 
ulate the observing process and produce simulated observa- 
tions that match the resolution and expected instrumental 
characteristic of the chosen camera. As described in Sec- 
tion [H we have used the SPIRE simulator to produce simu- 
lated observations for the 250-500 /xm data. 

Figure 8 shows the model IRDCs before (blue colour 
table) and after (red colour table) they have been processed 
into simulated observations. The columns are labelled at the 
top with the wavelength; the rows are labelled on the right- 
hand side with the description of the cloud. The maps show 
a radius of 1.1 parsec around the centre of each cloud. The 
pixel size of the model cloud maps at all wavelengths is 0.019 
pc (1.9 arcsec) for the spherical and flattened cloud, and 
0.028 (2.9 arcsec) for the turbulent cloud. The pixel sizes of 
the 250, 350, and 500 /im maps are approximately a third 
of the SPIRE beam FWHMs, corresponding to 0.058 pc (6 
arcsec), 0.097 pc (10 arcsec), and 0.116 pc (12 arcsec) at 
250, 350, and 500 /im respectively. The 8- and 15-/im bands 
are not shown in Figure 8, as they simply show up as dark 
clouds against a bright background. 

IRDCs have been identified in large numbers from infra- 
red surveys (Simon et al. 2006; Peretto & Fuller 2009) as 
higher-mass and more distant analogues to the optically 
dark clouds of gas and dust seen in silhouette against a 
bright background in V-band surveys (e.g. Barnard 1927; 
Lynds 1962). Starting at the shortest wavelength in Fig- 
ure 8 (70 /xm) this extinction can be seen in the spherical 
and flattened model clouds as a central dark patch. The 
size of this patch decreases with increasing wavelength un- 
til 500 /im where it has either closed up entirely or can no 
longer be resolved. We note that the size of the emitting 
region at 500 /xm is approximately the same size as the dark 
region at 70 /xm. 

The substructure of the turbulent cloud which is evident 
in the model images is hardly visible in the simulated images 

^ PACS has been developed by a consortium of institutes 
led by MPE (Germany) and including UVIE (Austria); KU 
Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA 
(Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); 
lAC (Spain). This development has been supported by the 
funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), 
CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and 
CICYT/MCYT (Spain). 



6 D. Stamatellos et al 

lOfim 160/zm 250/xm 550fj,m 500/xm 




Figure 8. A grid of maps of the three model IRDCs shown by wavelength. The wavelength is listed at the top of each column. The 
first, third, and fifth rows show the intensity maps from the radiative transfer modelling for the spherical, flattened, and turbulent clouds 
respectively. The second, fourth, and sixth rows show the map from the row above (250-500 /im only) after it has been processed by the 
SPIRE telescope simulator. The RT models are shown with a blue colour table, the simulated observations are shown with a red colour 
table. The individual maps are linearly scaled between zero intensity (black) and maximum on the map (white). 




Figure 9. Radial profiles of the three model IRDCs at wavelengths 8 — 850 |xm (top row). Flux from the simulated observations (solid 
lines and error bars) versus the RT model flux (dashed lines) at the SPIRE bands (middle row). SED fitted temperature using just the 
SPIRE bands (blue curve) and the SPIRE bands plus 160 fim (purple curve). The equivalent fit using the three SPIRE simulated maps 
is shown by the error bars. The dots correspond to the actual dust temperature computed by the radiative transfer models (bottom row). 



8 D. Stamatellos et al. 



Table 2. The geometry of the plane-of-t he-sky intensity distri- 
bution towards each of the model clouds. 



IRDC 


Observed FWHM (pc) 




250/im 


350 /im 


500 /im 


Spherical 


0.563 X 0.565 


0.474 X 0.474 


0.262 X 0.263 


Flattened 


0.565 X 0.368 


0.509 X 0.331 


0.341 X 0.218 


Turbulent 






1.31 X 1.31 



due to the angular resolution of the simulations being larger 
than the length scale of the substructure at the assumed 
distance of the IRDC (2 kpc). 

Table Elhsts the FWHM of 2D Gaussians fitted to each 
of the 250-500 /xm model clouds. These show quantitatively 
that the FWHM of the emitting region also decreases with 
increasing wavelength. This can be ascribed to the centre-to- 
edge temperature gradient. Cold dust located towards the 
centre of the cloud will preferentially radiate at longer wave- 
lengths, whereas warmer dust located towards the edge of 
the cloud will preferentially radiate at shorter wavelengths. 
The wavelength dependent size of the cloud is also visible in 
the simulated observations. However, the central extinction 
hole is only resolved for the spherical cloud at 250 /xm. 

5.2 Radial Profiles 

Figure 9 shows radial profiles of each of the IRDCs shown in 
Figure 8. Each column shows a single model cloud as anno- 
tated in the upper plot (spherical, flattened, or turbulent) 
while each row of panels represents a particular property 
or set of parameters. No a priori knowledge (e.g. about the 
density profile and the dust temperature) is assumed whilst 
generating the radial profiles. 

We compute an equivalent semi-major axis ai for each 
pixel i with coordinates {xi,yi) by assuming that it lies on 
an ellipse with the same origin, position angle, and axis ratio 
rsso as the equivalent 350 /xm cloud (i.e., the simulated ob- 
servations use the axis ratio of the 350 /xm simulated cloud) . 
The spherical and turbulent models have axis ratios of 1.0 
and the flattened model cloud has an axis ratio of ^ 1.5; it is 
assumed that the clouds' axes have a position angle of zero. 
The equivalent semi- major axis ai is given by 

= (Xcen - Xi)^ + (r350 X (^/cen " Vi))^ (4) 

where ( ^^cenj^/cen) are the coordinates of the centre of the 
cloud. Once ai is known, we bin the pixels into bins of width 
Aa and compute the mean and standard deviation of the 
intensities of the pixels in each bin. For the model clouds 
we use a bin width equal to the pixel size. This is 2.9 arcsec 
(0.029 pc). 

The top row of Figure 9 shows the distribution of mean 
intensity versus equivalent semi-major axis for each of the 
eight model bands. The flux density has been normalised for 
clarity, but no minimum flux has been subtracted (i.e. zero 
on the relative scale equates to zero on the absolute scale). 
The colour key for each band is shown beneath Figure 9. 
All three clouds show the characteristic IRDC extinction 
holes at PACS equivalent wavelengths (70 and 160 /urn). The 
spherical and flattened clouds are completely extinguished 
at 0.2-0.3 pc, but the turbulent cloud is extinguished across 
its entire radius (due to its higher mass). 



The middle row of Figure 9 shows the profiles of the 
three SPIRE bands. The model profile, from the upper row, 
is reproduced as a dashed line. Plotted over this is the profile 
of the simulated observation. The process of observing the 
clouds causes flux to be averaged out due to the resolution 
limit of the telescope. This moves flux towards regions of 
lower flux, i.e. outward or towards central holes, and is most 
pronounced for bands where there is a significant central dip 
in the model profiles. 

The lower panel of Figure 9 shows the results of fit- 
ting an optically-thin 13 — 1.85 (Ossenkopf & Henning 1994) 
greybody to the radial spectral energy distribution of each 
cloud. The blue and purple curves show the best-fit tem- 
perature for a greybody fitted to the 250-500 /xm and 160- 
500 /im model cloud profiles respectively (i.e. using the re- 
sults of the radiative transfer model). The error bars show 
the greybody fit to the 250-500 /xm simulated observation 
profiles (i.e. taking into account instrumental effects). These 
are interpolated to the lowest common resolution (500 /xm) 
before fitting the greybody. 

These figures show that temperatures fitted to the ob- 
served fluxes match the temperatures that would have been 
fitted to the model IRDCs in the absence of instrumental 
effects. However the fitted temperatures are higher than the 
actual temperatures (i.e. the ones calculated with the ra- 
diative transfer model) in the centre of the cloud and lower 
than the actual temperatures at the outer parts of the cloud. 
This is because the dust temperature from the simulated 
observations corresponds to an averaged temperature in the 
observed column, just as in real observations. 

This demonstrates that for accurately determining the 
dust temperatures in IRDCs, detailed radiative transfer 
modelling is needed along with mult i- wavelength, spatially 
resolved observations. Accurate dust temperatures are im- 
portant for determining the masses of these clouds. Previous 
studies have shown that overestimating cloud dust temper- 
atures even by a few degrees may lead to underestimating 
cloud masses by a factor of 2 — 3 (Stamatellos et al. 2007). 
This has important implications for the inferred cloud sta- 
bility and for the derived mass function of IRDCs. 

The methodology to determine the actual properties of 
an IRDC is to fit the observed fluxes at all available wave- 
bands, by varying the assumed density structure of the ra- 
diative transfer model. The density structure can be spher- 
ically symmetric or axisymmetric (e.g. flattened core), de- 
pending on each specific cloud. 3-dimensional modelling is 
also possible but the determined density structure cannot 
uniquely recovered (e.g. Steinacker et al. 2003, 2005). Once 
the observed fluxes are fitted, the density profile of the cloud 
is determined, and the actual dust temperatures are calcu- 
lated through the radiative transfer model. 



6 G29.55+00.18, AN IRDC FROM THE HI-GAL 
SURVEY 

G29. 55+00. 18 (Simon et al. 2006) is an infrared-dark 
cloud, which has been observed by Herschel/P ACS and 
Herschel/ SPIKE as part of the Hi-GAL survey, and by 
Spitzer/ GLIMPSEj. The cloud is seen in absorption against 
the bright background at 8 \im and 70 ^im, and in emission 
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Table 3. The physical properties of the infrared-dark cloud 
G29. 55+00. 18. i^cloud^ radius, ro : flattening radius, ric (H2) : 
central density, Tq_qqO /tq_qo: ratio of the visual optical depths, 
/isrf: ISRF factor (see text), Tsed: temperature calculated from 
a single-temperature, grey-body fltting of the SED, Tmodel^ tem- 
perature from the radiative transfer model, MgED- mass calcu- 
lated from Hildebrand (1983) and Tsed, ^modei: mass calculated 
from the radiative transfer model. The distance is taken from 
Heyer et al. (2009). 



G29.55+00.18 



R.A. (2000) 


18:44:37.00 


Dec. (2000) 


-02:55:07 


Distance (kpc) 


4.8 


^cioud (pc) 


1.7 


ro (pc) 


0.2 


nc(H2) (cm-3) 


2.8 X 10^ 




1.55 


/iSRF 


2.5 


TSED (K) 


16 


^model (K) 


10 - 21 


MsED (Mq) 


530 


Mmodel (M0) 


520 



at longer wavelengths (Fig. flO)) . The SED of the cloud is 
shown in Fig. 1121 

Guided by the image of the cloud at 500 |J,m, we model 
this cloud using a 2D density profile (Eq. [2]) with a ratio 
of optical depths T^_^^o/r^_^o = 1.55. The input parame- 
ters of the model and the derived physical properties of the 
cloud are listed in Table O The mass of the cloud is con- 
strained by the longer wavelength data, which are approxi- 
mately column density profiles, as the effect of the temper- 
ature gradient is relatively weak. The temperature of the 
cloud is mainly constrained by the shorter wavelength data. 
To match these short wavelength data, the external radi- 
ation field is enhanced by a factor of /isrf = 2.5, when 
compared with the ISRF at the solar neighbourhood, which 
is consistent with the higher ambient radiation field in the 
Galactic plane. 

The model reproduces well the images of the cloud at 
different wavelengths fFig. [TT]) . At shorter wavelengths (8, 
70 i^m) the cloud is seen in absorption (note that the back- 
ground is not modelled, hence there is no actual correspon- 
dence with the observed background), and in emission at 
longer wavelengths (> 170 ^im). The fiux and the shape of 
the cloud are well matched (assuming a viewing angle of 
^obs = 40°). At 170 \im there is secondary peak that is not 
reproduced by the model. This is probably due to asym- 
metric heating as this secondary peak does not appear at 
longer wavelengths (cf. Nutter et al. 2009); such heating is 
not included in the current model. 

The observed SED is well fitted by the model (Fig. [H 
dotted line). The dust temperature inside the cloud drops 
from to 21 K at its edge to 10 K at its centre (Fig. [131 cf. 
Peretto et al. 2010). Additionally, there is a variance of the 
temperature with the polar angle 0, with the more dense 
region at the "equator" of the cloud being colder by up to 
^ 3 K, as compared with the corresponding equal-radius 
region at the "pole" of the cloud. The mass of the cloud is 
calculated to be 530 Mq. 

We also compare the radiative transfer model with a 




1 00 1 000 

A (/xm) 



Figure 12. SED of G29. 55+00. 18. The dotted hue corresponds 
to the SED calculated from the radiative transfer model, and the 
dashed line corresponds to the SED from a single-temperature, 
grey-body flt. 




0.5 1 1.5 

r (pc) 



Figure 13. Dust temperature proflle of G29. 55+00. 18 at two 
different directions in the cloud. The bottom line corresponds 
to the midplane of the flattened structure {0 = 90°), whereas 
the upper line corresponds to the direction perpendicular to the 
midplane (6* = 0°). 



simple grey-body single-temperature fitting of the SED. The 
SED is well reproduced using a temperature of 16 K fFig. 1121 
dashed line). The mass of cloud calculated using the above 
temperature and the fiux at 500 ^im, is 520 Mq, which is 
slightly lower but consistent with the mass obtained from the 
radiative transfer model, despite using an average, and not 
the actual, dust temperature. However, the inner dense re- 
gions of the cloud where star formation will occur are colder 
(by ^5 K) than the estimated average temperature. 

This example demonstrates the need for radiative trans- 
fer modelling for determining the temperature structure of 
IRDCs. Despite the fact that a simple grey-body single- 
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Figure 10. G29. 55+00. 18 observations at 6 wavebands (top row: 8 |J.m Spitzer/GIAMFS^^ 70 |xm and 170 |J.m PACS; bottom row: 
250 |J.m, 350 |J.m and 500 |J.m SPIRE). The cloud is seen in absorption against the bright background at 8 |xm and 70 |xm, and in emission 
at longer wavelengths. 




Figure 11. G29. 55+00. 18 simulated observations at 6 wavebands (same as in Fig. llOj) . Contours are taken from the observed data at 
the corresponding wavelength. The model reproduces the general characteristics of the appearance the core (flux and shape). 
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temperature model provides a good fit to the SED, this 
model cannot accurately determine the temperature struc- 
ture in the cloud and overestimates the temperature of the 
interesting, in terms of potential for star formation, inner 
dense regions of the cloud. 



7 CONCLUSIONS 

We have demonstrated the use of the 3D, multi- wavelength 
Monte Carlo code phaethon, to model the transfer of ra- 
diation in IRDCs, that are externally illuminated by the 
interstellar radiation field. The cores of these clouds are be- 
lieved to be where high-mass stars form, i.e. they are the 
high- mass equivalent of prestellar clouds. We have presented 
three widely different models, in which we varied the mass, 
density, radius, morphology and internal velocity field of the 
cloud. We have shown the predicted output of the model at 
the wavebands of Herschel and Spitzer. We also passed the 
model output through the SPIRE simulator to produce sim- 
ulated observations of these IRDCs. These were then anal- 
ysed as if they were real observations. Subsequently, the re- 
sults of this analysis were compared with the results of the 
radiative transfer modelling. 

Our study highlights the need for detailed radiative 
transfer modelling when using multi-wavelength observa- 
tions from Herschel to accurately determine the proper- 
ties of IRDCs. This method was applied for the study of 
G29. 55+00. 18, an IRDC from the Hi-GAL survey. Modelling 
of a larger sample of IRDCs found in the Hi-GAL survey will 
appear in a future publication (Wilcock et al., in prep). 
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